function [Km2_1,Km2_2,Km2_3,Km2_4] = lpcvd_Km2_sincos(L,mode,method)
Km2_1 = zeros(mode+1,mode+1);
Km2_2 = zeros(mode+1,mode+1);
Km2_3 = zeros(mode+1,mode+1);
Km2_4 = zeros(mode+1,mode+1);

dx = pi/L;
x = 0:dx:pi-dx;
y = 0:dx:pi-dx;
if nargin == 2
    method = 'dis';
end
if strcmp(method,'dis')
    for m = 0:mode
        for n = 0:mode
            for ix = 1:L
                for iy = 1:L
                    if ix == L
                        x_next = x(1);
                    else
                        x_next = x(ix+1);
                    end
                    Km2_1(m+1,n+1) = Km2_1(m+1,n+1)+(fun_phi1mn(m,n,x_next,y(iy),L)-fun_phi1mn(m,n,x(ix),y(iy),L))^2;
                    Km2_2(m+1,n+1) = Km2_2(m+1,n+1)+(fun_phi2mn(m,n,x_next,y(iy),L)-fun_phi2mn(m,n,x(ix),y(iy),L))^2;
                    Km2_3(m+1,n+1) = Km2_3(m+1,n+1)+(fun_phi3mn(m,n,x_next,y(iy),L)-fun_phi3mn(m,n,x(ix),y(iy),L))^2;
                    Km2_4(m+1,n+1) = Km2_4(m+1,n+1)+(fun_phi4mn(m,n,x_next,y(iy),L)-fun_phi4mn(m,n,x(ix),y(iy),L))^2;
                end
            end
            Km2_1(m+1,n+1) = Km2_1(m+1,n+1)/(pi^2);
            Km2_2(m+1,n+1) = Km2_2(m+1,n+1)/(pi^2);
            Km2_3(m+1,n+1) = Km2_3(m+1,n+1)/(pi^2);
            Km2_4(m+1,n+1) = Km2_4(m+1,n+1)/(pi^2);
        end
    end
elseif strcmp(method,'ana')
    for m =0:mode
        for n = 0:mode
            if m == 0 || 2*m==L || n==0 || 2*n == L
                Km2_1(m+1,n+1) = 0;
            else
                Km2_1(m+1,n+1) = 4/(pi^2)*(sin(pi*m/L))^2;
            end
            
            if m== 0
                Km2_2(m+1,n+1) = 0;
            else
                Km2_2(m+1,n+1) = 4/(pi^2)*(sin(pi*m/L))^2;
            end
            
            if m==0 || 2*m==L
                Km2_3(m+1,n+1) = 0;
            else
                Km2_3(m+1,n+1) = 4/(pi^2)*(sin(pi*m/L))^2;
            end
            
            if n==0 || 2*n==L || m == 0
                Km2_4(m+1,n+1) = 0.0;
            else
                Km2_4(m+1,n+1) = 4/(pi^2)*(sin(pi*m/L))^2;
            end
        end
    end
else
    error('Invalid method');
end

    function phi1 = fun_phi1mn(m,n,x,y,L)
        phi1 = (2/L)*sin(2*m*x)*(sin(2*n*y))';
    end

    function phi2 = fun_phi2mn(m,n,x,y,L)
        if (m== 0 || 2*m == L) && (n==0||2*n==L)
            coeff = 1/L;
        elseif m~=0 && 2*m~=L && n~=0 && 2*n~=L
            coeff = 2/L;            
        else
            coeff = sqrt(2)/L;
        end
        phi2 = coeff*cos(2*m*x)*(cos(2*n*y))';
    end

    function phi3 = fun_phi3mn(m,n,x,y,L)
        if n == 0 || 2*n == L
            coeff = sqrt(2)/L;
        else
            coeff = 2/L;
        end
        phi3 = coeff*sin(2*m*x)*(cos(2*n*y))';
    end

    function phi4 = fun_phi4mn(m,n,x,y,L)
        if m ==0 || 2*m == L
            coeff = sqrt(2)/L;
        else
            coeff = 2/L;
        end
        phi4 = coeff*cos(2*m*x)*(sin(2*n*y))';
    end

end